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VARIABLE VISCOSITY AND DENSITY BIOFILM SIMULATIONS USING AN 
IMMERSED BOUNDARY METHOD, PART II: EXPERIMENTAL VALIDATION AND 

THE HETEROGENEOUS RHEOLOGY-IBM 


JAY A. STOTSKY* * * § , JASON F. HAMMONDt, LEONID PAVLOVSKY*, ELIZABETH J. STEWART*, JOHN 
G. YOUNGER§, MICHAEL J. SOLOMON*, AND DAVID M. BORTZ*’ 

Abstract. The goal of this work is to develop a numerical simulation that accurately captures the biomechanical response 
of bacterial biofilms and their associated extracellular matrix (ECM). In this, the second of a two-part effort, the primary 
focus is on formally presenting the heterogeneous rheology Immersed Boundary Method (hrIBM) and validating our model 
against experimental results. With this extension of the Immersed Bounadry Method (IBM), we use the techniques originally 
developed in Part I, (Hammond et al. El) to treat the biofilm as a viscoelastic fluid possessing variable rheological properties 
anchored to a set of moving locations (i.e., the bacteria locations). We validate our modeling approach from Part I by comparing 
dynamic moduli and compliance moduli computed from our model to data from mechanical characterization experiments on 
Staphyloeoeeus epidermidis biofilms. The experimental setup is described in Pavlovsky et al. (2013) |22| in which biofilms 
are grown and tested in a parallel plate rheometer. Matlab code used to produce results in this paper will be available at 
https: / / github.com/MathBioCU/BiofilmSim. 
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1. Introduction. The goal of this work is to develop a numerical simulation method that accurately 
captures the biomechanical response of bacterial biofilms and their associated extracellular matrix (ECM). 
In this second paper, we show that the model and simulation method developed in part im, can be used to 
predict material properties of a biofilm and that the simulated results mimic experimentally measured results. 
The underlying mathematical technique is an adaptation of the Immersed Boundary Method (IBM) that 
takes into account the finite volume of bacteria, and variable material parameters found in biofilms whose 
variation is anchored to the positions of bacteria in a biofilm. We call this method the heterogeneous rheology 
Immersed Boundary Method (hrIBM). A key feature of our results is that the simulations are initialized with 
experimentally measured position data providing the locations of bacteria in live S. epidermidis biofilms. 
This removes ambiguity about how to represent the biofilm computationally. When using this data, the bulk 
physical properties estimated through simulation match experimental results. We also verify that when using 
different position data sets that possess similar spatial statistics, the physical properties of the biofilm do not 
change significantly. We also provide quantitative results on the periodic rotation of suspended aggregates 
of bacteria in shear flow. 

In recent years, much work has been done to develop detailed mathematical models that capture the 
biomechanical response of bacterial biofilms to physical changes [D El 0 [la [IS]. In general, the physical 
properties governing the growth, attachment, and detachment of a biofilm are dependent on the ECM, a 
viscous mixture of polysaccharides and other biological products excreted by bacteria in the biofilm. The 
focus of this work is on accurately simulating the biomechanical response of a biofilm and its associated 
ECM due to applied shear stress and shear strain. 

In Section]^ we provide a brief review of the classical Immersed Boundary Method (IBM), a well known 
computational technique used for the simulation of coupled fluid-structure interactions. Additionally, we 
discuss some other IBM based biofilm models, and explain the adaptations of the IBM that lead to the 
hrIBM. In Sectiona description of the numerical properties and, results from numerical tests showing that 
the model is convergent are provided. In Sectionmethodologies for computing relevant material properties 
from the model are discussed, and the dynamic moduli and compliance moduli estimated by the model are 
compared to experimental data from biofilms grown in a bioreactor. We observe that these properties do not 
significantly vary when several different experimental coordinate data sets with similar spatial statistics are 
used. We also compare results of tumbling of bacteria aggregates suspended in shear flow against theoretical 
results provided in Blaser et al. [3]. In Section ([^, we discuss future research direction and limitations. 
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The ability to calculate bulk material properties of a biofilm while directly incorporating the microscale 
rheology and connectivity of the biofilm is the primary contribution of this article. This development shows 
that IBM-based models which connect fine scale features such as models that describe viscoelastic connections 
between bacteria, to fluid dynamical models can provide physically accurate results. From our results, we see 
that the hrIBM model accurately captures the elastic component of the biomechanical response of biofilms 
to applied stress and strain, and matches experimental trends observed in the viscous response. 

To our knowledge, this work is the first to use a model that accounts for both the heterogeneous rheolog¬ 
ical properties and the inter-bacterial connectivity to compute material properties of a biofilm. Code used to 
produce the results obtained in this paper will be available at https://github.com/MathBioCU/BiofilmSim. 

2. The Biofilm Model. In this section, we discuss some previous biofilm models and explain the 
alterations of the classical IBM that lead to the hrIBM. In Section we introduce our biofilm model. In 
our model, we couple a spring model of the inter-bacteria links in the biofilm with fluid motion through 
the biofilm to treat the biofilm as a multicomponent viscoelastic material. On the level of our simulations, 
both the fluid-structure interactions of the bacteria and the surrounding fluid, and the interconnectedness 
of bacteria in the biofilm play a major role. 


2.1. Previous IBM Based Biofilm Models. In recent years, a number of different approaches to 
IBM-based biological material models have been developed. One such biofilm model can be found in Luo 
et al. m- In their model, they couple an immersed viscoelastic structure to the fluid flow in an immersed 
boundary type formulation. However, the fluid equations are solved separately from the equations governing 
the motion of the immersed viscoelastic solid and then coupled together at a physical interface. Our model 
builds upon this work by eliminating the need for an explicit interface since biofilms frequently do not have 
well defined fluid-structure interfaces. Another approach to capturing the viscoelastic nature of biofilms with 
an the immersed boundary method is through the choice of viscoelastic model used for the links between 
bacteria. The choice of model can be used to affect the value of the external force density, f, in the Navier- 
Stokes equations, (2.1). This type of strategy was used first by Bottino in [4] to model general viscoelastic 


connections in actin cytoskeleton of ameboid cells and also by Dillon and Zhuo in [28] to model sperm 
motility. 

IBM-based models can be found in Alpkvist and Klapper, [2] and in Dan Vo et al. [9|. In these models, 
an IBM is used directly to couple the forces between connected bacteria with fluid motion. Additionally, 
some validation results are performed to show that properties such as the recovery and relaxation times of a 
biofilm can be modeled with such a method. Our model also builds on these by including spatially variable 
rheological properties which are an important structural feature of biofilms. Additionally, we note the recent 
work of Sundarson et al. [26] in which an IBM model is used to model detachment of biofilms. 

We would also like to point out that detailed explanations of the IBM can be found in usi m im ng, 
and additional IBM-based biofilm models can be found in [mails]. 

2.2. The Biofilm Model. The model we use is comprised of two sets of equations; those that model 
fluid flow through the biofilm, and those that model motions and forces experienced by each bacteria cell in 
the biofilm. These equations are listed in 

( 2 . 1 ) 

( 2 . 2 ) 

(2.3) 

(2.4) 

(2.5) 

(2.6) 

(2.7) 

( 2 . 8 ) 


p{x, t) (ut -I- (u ■ V) u) = -VP -I- V ■ /n(x, t) (Vu -1- (Vu)^) -f f(x, t) 

V-u = 0 

U(X(s, t), t) = f u(x, t) S{'X{s, t) — x) dx s = 1, 2, ...N 

Jn 

9X(s, t) _ . . . . 

— U(X(5, t), t) 

F(X(5, t),t)=.F(X(5, t),V) 

f (x, t) = ^ [ F(X(s, t), t) (5(X(s, t) - X, u) dX 

p{x, t) = po + min | J ui^pi, S{X{s, t) - x, lo ) dX, pbj 

p(x, t) = po + min < / LO^pb <^(X(s, t) — x, lo ) dX, pb 





Symbol 

Definition 

u, P, F 

Eulerian velocity, pressure and force density 

p{x, t), p{x, t) 

spatially (and temporally) varying density and 
viscosity 

coo 

Average volume associated with a cell and its 
surroundings, number of bacteria in the domain 

X(s, t), U(-, f), F(., t) 

Bacteria position, velocity, and force density, 
labelled by Lagrangian coordinate, s 

O 

o 

Density and viscosity of pure water 

ph, Pb 

Density and viscosity of the biofilm at the center 
of mass of each bacteria 

H; V) 

Function that determines the force associated 
with each bacteria based on a constitutive 
viscoelastic model, V 


The Dirac delta function 

(5(-, w) 

A smoothed approximation of the Dirac delta 
function. The 2nd argument, cj is a 
hydrodynamic parameter corresponding to the 
radius of a bacterium. 


Table 2.1: List of terms in governing equations. We attempt to use standard notations when possible. 


The same set of equations was also used in Part I, m and are reproduced here for the convenience of the 
reader. The quantities appearing in these equations are listed in Table |2T| 

Since individual bacteria are not assumed to have infinitessimal volume at the scale of our simulations, 
the Lagrangian quantities; X, F, a nd U correspond to measurements taken at the center of mass of each 
bacterium. As described in Section 2.3, the second argument, cj, of the smoothed Dirac S function, <5(-, uj) 
determines a region of support for the smoothed S function . The choice of ^(•, uj ) govern how the mass 
density, viscosity, and force density vary around each bacterium. 

Additionally, since in each simulation, there is a fixed number, N of bacteria in the computational 
domain which is independent of the grid spacing, h, we can write the integrals in equations (2.6)-(2.8) as 
summations of the form 


f(x, t) = 


1 

^^F(X(s, t),t)5{Xis, t)-x, w) 

^ S = 1 


N 


p{x, t) = po + min < '^L0^pb5{X{s, t) - x, w), 


Pb 


< s = l 


f ^ 

p{x, t) = po + min < y] (5(X(s, t) - x, to), pt, 

U=i 

These summations are slightly different than those used on part I. 

By using the IBM as a basis for our biofilm model, we are able to avoid treating the biofilm as a two phase 
fluid with a distinct bulk fluid region and a distinct biofilm region. Instead, the use of variable rheological 
properties over the entire domain couples the biofilm and bulk fluid motions as a single viscoelastic material. 

2.3. The Heterogeneous Rheology Immersed Boundary Method. We will now describe our 
reasoning behind equations ( |2.1| )-( [2^ . In our model, we extend the IBM to account for the fixed, finite size 
of bacteria and allow for variable physical properties that are anchored to a moving Lagrangian mesh (i.e. 
the bacteria positions). We denote this approach the heterogeneous rheology IBM (hrIBM). 
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Lagrangian Variables 



Eulerian Variables 


Figure 2.1: The coupling between the Eulerian and Lagrangian variables in the hrIBM is shown here. The 
Eulerian and Lagrangian variables are coupled by the couiputation of U from u, and the computation of f , 
/i, and p from F and X. The IBM is a widely applicable method in part because it allows for a great variety 
of fluid solvers and solid structural models to be coupled through 5 function transfer identities. 


The original IBM was first developed as a means of solving fluid-structure interaction problems in 
cardiology and is applicable to problems with moving, irregularly shaped boundaries [531 US]- With the 
IBM, the fluid velocity fields and pressure are usually solved for on a fixed, Eulerian grid and the movement 
of the boundaries due to fluid motion is tracked by a moving Lagrangian mesh. As material boundaries are 
deformed, a constitutive model is used to determine the force density exerted by the boundary on the fluid 
around each Lagrangian point. The Lagrangian force density field is then transferred to an Eulerian force 
density field through the use of a discrete approximation of the following identity, 

(2.9) f (x, t)= [ F(X(q, t)) VX(q, t) - x) dq 

Jn 


where d(-) is the Dirac S function. In our biofilm model, the transfer of quantities from the Lagrangian to 
the Eulerian grid is done with a smoothed S function that differs from the standard choices used in most 
IBM literature (see ( |2.2[ )). The effect of the Eulerian force term on the velocity and pressure fields is then 
found by solving the Navier-Stokes equations, with appropriate boundary conditions. 

In the original IBM, after discretizing, the integration in (2.9) is carried out by computing a sum of the 
form 


( 2 . 10 ) 


N 

f(xj, tj) = y]F(Xfc, tj)6{Xk -Xj, h) h^, 
k 


where the Eulerian and Lagrangian forces are evaluated at the Eulerian and Lagrangian grid points respec¬ 
tively, and <5(-, h) is a discrete approximation of the Dirac delta function that has compact support related 
to the grid spacing parameter h. With the IBM, the discrete approximation is chosen such that as h ^ 0, 
d(r, h) d(r). This makes sense for fluid structure interactions involving fluid-solid boundaries that have 
infinitessimal thickness, and thus zero volume. In biofilm modeling, each Lagrangian point corresponds to 
















the center of mass of a bacterium which has finite dimensions. Therefore, we use a smoothed version of the 
standard discrete S function that has a fixed region of support, independent of the grid spacing, which is 
governed by a radial parameter, cj. 

In our model, we use a smoothed discrete Dirac S approximation of the form, 

(2.11) j(x,„) = ±,^(£),f(i)0(£) 

with 0(r) as defined in [24] by 

{ I (5 - 2 |r| - ) l<kl<2 

I ^3 — 2|r| + yi + 4|r| — 4|rP^ 0 < |r| < 1 

0 |r|>2 

This is chosen because it most closely satisfies the unity and first-moment conditions described below 
for the values of cj we use. If cj = h, the standard discrete S functions seen in IBM literature is obtained. For 
this work, we assume that the bacteria are spherical and thus cj is understood as a hydrodynamic radius. 
We also note that extensions to this formalism will allow for the treatment of nonspherical bacteria. Thus, 
cj may be though of more generally as a shape parameter. 

With <5(x, h), the unity condition, 

(2.12) y] (5(x-X,= 1, VX, 

xeGh 


and first-moment condition, 


(2.13) 




(x - X) (5(x - X, h)h^ = 0, VX, 


are both satisfied. With a grid-independent choice for cj, these properties are only satisfied approxim ately . 
However, we do see that in the limit as h ^ 0, greater than O(h^) convergence in ^(r, uj) to equations (2.12) 


and (2.13) is observed. 


Highly heterogeneous viscosity and moderately heterogeneous density are common characteristics of 
biofilms. Although IB methods with variable density have existed for some time (see m), the incorporation 
of spatially variable viscosity in the IBM is an area that has yet to be well developed. We do, however note 
the recent publications by Fai et al. mm in which an IBM capable of solving problems with variable 
viscosity and density is used to model the motion of red blood cells flowing in capillaries. When modeling 
red blood cells, the viscosity exhibits a “jump” discontinuity between the blood plasma and the intracellular, 
hemoglobin-containing fluid of a red blood cell. Thus, their model is designed to capture the dynamics of 
two interacting fluids with different rheological properties separated by a deformable membrane. In our case, 
there do not exist well defined boundaries and thus, (!(•, uj) is adjusted to reflect this. 

In biofilms, the spatial variance of material properties is localized around the position of each bacterium, 
while in fluid far away from any bacteria, the physical properties are those of the bulk fluid. This localization 
of the variation in material properties allows the spatial variation in density and viscosity to be found by 
using a smoothed d-function integration similar to that used to compute the Eulerian force field. We define 
an effective viscosity, and an effective density, and assume that at the center of mass of each bacteria, 
the viscosity and density are p(X^, t) = p^, and t) = p^. Defining po fo be the viscosity of the bulk 

fluid, in this case water, the viscosity at any Eulerian grid point can be calculated as: 


(2.14) 


//(x) = min Mo + / w^(M 6 “ Mo) “ X(s), uj) dX, ^ib 

Jq 


A similar formula exists for the spatial variation of density. A summation sign is used instead of an integral 
since the number of bacteria, N, is fixed and independent of the mean Lagrangian mesh spacing do. 

In this model, we indirectly take into account the fluid volume displacement caused by the presence of the 
bacteria. We treat the localized high viscosity around each bacteria as an effective viscosity that accounts for 
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Bacteria Positions 




Bacteria Positions 


Isosurface, = 125//o 



(c) 


Figure 2.2: a) Shows the 3D locations of bacteria from experimental biofilm data, b) Each line represents a 
viscoelastic connection between two bacteria. Bacteria connected if they are within 1.62/im of each other, 
c) A viscosity isosurface of the same biofilm. The maximum viscosity is 250/io where /io is the viscosity of 
water. The isosurface is the surface defined by /i(x) = 125/io- 


both the displaced fluid volume and the increased viscosity near the bacteria surface m Extensions based 
on changing our choice for S(x, cj) could possibly allow for a more precise computation volume displacement 
into the model. The bacteria 5. epidermidis is known to have a diameter of approximately 0.5 — 1.0 /rm [13] 
, thus we choose uj such that the viscosity halo around each bacterium is a little greater than 1 /rm in our 
simulations. 


3. Numerical Methods. The numerical methods we use are based on those originally discussed in 
Hammond et al. m- We summarize them here for convenience and also provide convergence results. To 
approximate solutions to equations -( |2.8[ ), we use a projection method similar to that used in Zhu et 
al.[17|. The solution scheme uses an implicit Euler solver to update an intermediate velocity profile at each 
time step and is expected to be 0{At) convergent. To discretize the domain, we use a uniform finite difference 
discretization with equal spacings in the x, y, and 2 : directions. The spatial derivatives are approximated 
with 2nd order, centered finite differences. 


3.1. Numerical Algorithm. At each time step, the following quantities must be updated: u, U, P, 
F, f, X, /i, and p. To improve numerical accuracy, we nondimensionalize the problem with the following 
choices: 


U - P 'X f. P ^ t 

u= — P = ^ f=J- M = — P = — x= - t = — 

uo Po JO Po Po L to 

and also introduce the following nondimensional parameters: 


Re=P^ St=^ Cl 

Po toUo 


Po 

Poul 


C 2 


foL 

Poul' 


As is standard terminology. Re is the Reynold’s number, St the Strouhal number, and Ci and C 2 are 
additional constants. Additionally, we define dg to be the average Lagrangian volume element as described 
in Part I US]- Eor convenience, we will now assume that all quantities are nondimensional unless otherwise 
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Quantity 

Value 

Po 

IPa 

1^0 

1-10-^ Pa-s 

Po 

998 kg/nn? 

L 

10“^ m 

Re 

0(10-^) 

St 

0(10-'') 

Cl 

O(IO^) 

C2 

0(1) 

to 

1 s 

fo 

INjuP 

Uo 

0(10“^) (varies) 

do 

0.159 i 

Pb 

0.12po 

Pb 

250/io 

J- max 

1.3223 ■ 10“*^ 

Connection Distance 

0.162 P 

Radial Parameter, uj 

0.033 ■ L 


Table 3.1: Values of Physical Parameters and Nondimensional constants used in simulations 


stated. The values of the constants we use are listed in Table 13.11 and the motivation for these values is 
discussed in Part I. 

As is standard practice in IBM algorithms, we uncouple the Eulerian variable updates and Lagrangian 
variable updates for computational reasons. At each time step, we use a projection-based solver to solve 
the Navier-Stokes equation for u and P. We define Gh a discrete gradient operator, and Dh a discrete 
divergence operator, and use the following projection method to obtain u and P: 

1. Solve for u* 

(^St ^ 


2. Solve for P^’^) 


3. Compute 


Re 


D, 


(Gft(u*) + (Gft(u*))^)l 





fSt \ Uh{u*) 
\CJ At 


u^^'> = u* 



At 




In steps 1 and 2, full multigrid solvers and multigrid preconditioned conjugate gradient solvers are used to 
find u* and . After obtaining the updated velocity and pressure, the Lagrangian velocity and position 
updates follow. 


= Y, h) 

h^Qh 




At 


u("). 
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(a) (b) 

Figure 3.1: a) A viscosity isosurface is shown for a small section of a biofilm used in simulation. The inner 
isosurface is /i = 125/io and the outer transparent isosurface is at /i = 50/io- Slices of the ||u|| velocity field 
are shown as well, b) The eyz component of the strain rate is plotted on the /i = 125/io viscosity isosurface. 
Additionally, the strain rate and contours of viscosity are shown in slice planes. For a single phase fluid, 
Newton’s viscosity law is a = jie. Although biofilms are not Newtonian fluids, we still see that in areas of 
low viscosity, higher strain rates are found and in areas of higher viscosity lower strain rates occur. 


Next the Lagrangian force density is computed based on the new positions, as Finally, 

the Eulerian fields, q = {f , /r, p} are computed using discrete^ function interpolation to the Eulerian grid 
through equations of the form. 


q(") = w). 

x(’")er: 

In the simulations we conduct, the primary direction of fluid flow is in the ^ direction. The height is 
governed by the y coordinate and width by the x coordinate. 

3.2. Numerical Verification and Convergence Properties. In the first numerical verification re¬ 
sult, we verify the accuracy of the numerical projection method solver with no biofilm present by comparing 
the numerical solution with an analytical solution. Since there is no immersed structure, this is a test of the 
fluid solver alone, and not the IBM method. Eor this test, the domain, Vi is chosen to be a rectangular solid 
that is periodic in the x and ^ directions. Erom [3, the following boundary conditions for ^ = 0 and yL^ 

dP 

(3.1) u|o = 0 u|^^= (0, 0, sinz/t) 

provide us with an analytic solution. 


(3.2) 


Uz{y, t) 


sinh /c ^(1 + i) . / / sinhkyil i) \\ 

sinhkyL{l + i) y^(^i + i) j j 



The values of P, and Uy are exactly zero in this case. The values of p and p are set to 998kg/m^and 
1 Pa • 5 respectively and are homogenous across the domain since no analytic solutions with variable density 
and viscosity and the boundary conditions given above are known to the authors. Convergence tests were 














Frequency 

Time 

Space 

Error — u oo 

IHz 

1.006 

1.801 

4.910 • 10-i'J 

100 Hz 

1.221 

2.002 

2.223 ■ 10“*^ 


Table 3.2: Spatial convergence tests were carried out with grid spacings, h, set to 1/32, 1/64, and 1/128 
and a time step of At = 1/500/z/. Temporal convergence tests were done with vAt set to 1/125, 1/250, 
and 1/500 and Ax = 1/64. Error is computed at t = (0.2/z^)s. Convergence factors are computed as 

\\u{At/ 2 )-u{At )\\2 time and by p{h) = log where is an interpolation 


p{At) = log - 


|n(At/4)-n(At/2)||, - -5 | |n( V4)-//;,^n(M/2| ^ 

operator taking functions from a grid with spacing h to a, grid with spacing h/2. 


Frequency 

Velocity, u 

Position, X 

Time 

Space 

Time 

Space 

49.91Hz 

0.983 

1.105 

1.022 

0.952 

4.991 Hz 

0.991 

0.910 

1.007 

1.054 


Table 3.3: Convergence factors of hrIBM with biofilm. For spatial convergence, h was set to 1/32, 1/64, 
and 1/128 with a time step of uAt = 1/500. To measure the temporal convergence factors, uAt was set to 
1/250, 1/500, and 1/1000. In both cases, the boundary conditions described in Section were used. 


conducted with frequencies u = 1 Hz and u = 100 Hz. In Table [3i^ the absolute error, temporal convergence 
factors, and spatial convergence factors are listed. These quantities are calculated as described in Part I eg. 

Additionally, with the same boundary conditions as above, we tested the convergence rates for simulations 
with a biofilm that possesses variable density and viscosity. Temporal and spatial convergence factors are 
shown in Table |3.3[ More detailed numerical convergence results for this model with different boundary 
conditions are shown in m In Table [33| temporal convergence factors for the same fluid conditions and 


domain as the analytical solution are listed. The pressure convergence rate is not shown here since pressure 
variation only varies by about 0(10“^) and thus is around the same order as the numerical errors observed 
in the finite difference approximations used. 


4. Experimental Validation Results. The material characterization of bacterial biofilms is a difficult 
experimental task. It is usually not possible to grow biofilms large enough for use in standard testing devices 
and, attempts to move a biofilm from the environment it was grown in to a testing apparatus may alter its 
structure [22]. In Pavlovsky et al. (2013) [22], a promising experimental method of testing material properties 
of biofilms was developed. In the experimental setup, a biofilm is grown in a parallel plate rheometer. As 
the biofilm grows, it adheres to both the top and bottom plate of the device. The top plate can then be 
rotated or repositioned vertically and the stress and strain induced in the biofilm can be monitered. These 
measurements can then be used to infer material properties of the biofilm. Using the hrIBM model, we set 
up a simulation to reproduce experiments described in Pavlovsky et al (2013). 

In order to reproduce the biofilm in simulation, 3D position data sets obtained by high resolution 
microscopy of live biofilms are used to initialize the positions of bacteria in the computational domain. The 
experimental setup used to obtain these data sets are described in Pavlovsky et al. (2015) [21] and Stewart 
et al. [25]. Although the biofilm position data sets that we use, which were obtained from the experiments 
described in [21], are not the ones grown and tested in the bioreactor, they are from biofilms grown under 
similar physical and nutrient availability conditions. A key result seen from our simulations is that the 
material properties computed by our model of the different data sets are similar to each other. This indicates 
that the material properties obtained through simulation must depend on larger scale structural properties of 
the biofilm and may be treated as bulk properties of the biofilm. For validation we compared bulk properties 
measured by our model to experimental results. The methods used to compute these quantities are discussed 
in the next subsections. 

In Pavlovsky et al. (2013) [22], small amplitude rheometry (SAR) is used to characterize the viscoelastic 
behavior of S. epidermidis biofilms. In SAR experiments, the upper plate of the rheometer is rotated to 
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induce a sinusoidal shear deformation such that the average strain amplitude at the top of the biofilm is 
a fixed value and the corresponding stress is measured. The strain amplitude was set to 0.13 at the outer 
radius of the rheometer since this strain amplitude is found to be in a regime of primarily linear and elastic 
mechanical behavior [H]. Using 4.7, the dynamic moduli are computed at a number of different frequencies 
of oscillation. 

Creep compliance testing is another characterization technique used in Pavlovsky et al. (2013). In a 
creep compliance test, a constant shear stress is applied to the biofilm through the top plate of the rheometer. 
This induces a time dependent strain which can be measured. 

With the hrIBM model, we assume that for a small rectangular sample of the biofilm that is not near the 
rotational center and, does not border the outer boundary of the disc, the effects of cylindrical geometry are 
negligible and the rotational motion can be approximated as linear shear. This assumption greatly reduces 
the computational expense of simulating the biofilm and simplifies the discretization of the computational 
domain. This approximation is valid since the stresses due to angular momentum are much less than those 
due to the shearing motion of the plates. The size of the bioreactor used experimentally is 40 mm in diameter 
and approximately 250/im in height, whereas, the computational domain is only 9 — 18/im in width and 
length and 18 — 27 /im in height. 

From Christensen [S], the shear stress and strain, cjoz and eoz^ of an isotropic viscoelastic cylinder 
undergoing small angle torsion are proportional to r, the radial coordinate. Thus, if we choose to simulate 
some subset of the bioreactor that is 30 {im in width (this is larger than in simulations we conduct) that 
is near, but not touching the outer edge of the cylinder, at a radius of 15 mm from the center, the ratio of 
shear strain and shear stress exerted at the inner and outer boundaries is approximately (15mm)/(15mm + 
30 jam) = 0.998. Additionally, we see that in this case, the shear strain and shear stress are not functions 
of the angular coordinate, thus approximating the slightly curved domain as a rectangular solid should not 
alter the physics of the problem. Of course, biofilms have far more complicated material properties, however, 
at the scale of our simulations, this result indicates that rectangular geometry and linear shear produces an 
accurate approximation of the motion of the biofilm.. 


4.1. Computation of Rheological Properties. In order to compute the desired dynamic moduli, 
and compliance modulus results, the stress and strain experienced by the biofilm during simulation must 
be computed. The stress cr, is decomposed into a sum of stress due to the fluid motion, and stress 
due to the straining of inter-bacteria connections within the biofilm cr^. The total stress can then be found 
as tr = Although each component of stress is computed separately during simulations, distinct 

simulations cannot be used to individually test cr^ and since they are coupled. 

In order to calucate the strain e, a set of tracer particles is tracked throughout the simulation. Spatial 
derivatives can then be calculated to obtain approximations of the strain. Additionally, since only small 
amplitude strains are observed, the linear relation, e = |(Vd + V^d) is an accurate approximation of the 
strain for a displacement vector d. The derivatives needed to compute the strain are taken with respect to 
the advected material coordinates. 

Viscoelastic materials are often characterized through their time dependent stress response to strain or 
their time dependent strain response to stress. For a general viscoelastic material, given that the stress and 
strain are sufficiently smooth functions of time, constitutive relations between the stress and strain may be 
written in terms of a convolution with viscoelasticity tensors as: 


(4.1) 


crij(x, t)= Gijkliji-, t-r) 

J —OO 


dr 


efc;(x, T)dT 


(4.2) 


^) — / dijkiipc^ t t) 

J —OO 


dr 


(Tki{yi, r)dT 


where aij is the stress tensor, Cij is the strain tensor, and Gijki and Jijki are fourth order viscoelasticity 
tensors (see Christenson [8] , §1 for a derivation). In the literature, G is often called the relaxation modulus 
and J is called the compliance modulus. For linear, isotropic materials, the expression for Gijki simplifies 

to Gijki = I — G 2 {t)) 6ij6ki + ^Gi iSikSji -\-Sii6jk)^ where <5^^ is the Kronecker delta function and 
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Einstein summation notation is used. The two functions, Gi{t) and G 2 (t) correspond respectively to shear 
and dilatational stresses. Analogous expressions exist for the compliance tensor. Although the viscoelastic 
moduli are spatially heterogeneous, we believe that more meaningful results are obtained in the mean field, 
or spatially averaged, time dependent values for e, cr, G, and J. These quantities depend less on the exact 
configuration of bacteria in a biofilm and behave more like bulk material parameters that can be measured 
experimentally. Although the interconnected links used to model the connections between adjacent bacteria 
each individually introduce anistropy into the model, under the conditions of our simulations, the overall 
behavior of the biofilm is not highly anisotropic. 

4.1.1. Computation of Strain. Although a single phase Newtonian fluid will behave viscously (i.e., 
the stress only depends on the strain rate, not strain itself), in a biofilm the fluid component is influenced 
by the elastic components of the biofilm and thus the stress state in a biofilm depends directly on the strain 
(along with the strain rate). In order to compute the strain, the displacement field must be computed. The 
displacement, d of a particle located at xq at time, t = 0 in a material undergoing deformation can be found 
by solving the following ODE: 

d C 

(4.3) “ ^(^0’ 

dt Jq 

In the biofilm simulations, “tracer” particles with positions denoted by S(x, y, z)^ are initialized at heights 
Vl — Vl — 1 and — 7 — 2 h, near the top of the biofilm at t = 0 . At each time step, the positions of 
the tracers are updated using the same 6 function interpolation used to update the bacteria positions. With 
these tracers, the deformation of the biofilm can be tracked throughout the simulation. 

In the simulations, the eyz component of strain is needed at the upper boundary of the domain. Therefore, 
the tracers are initialized near the top of the domain in three vertically aligned layers. This is done to make 
the numerical approximation of derivatives of the form dd^jdy easier . With the initial arrangement of 
tracters in vertically aligned layers, the centered finite difference approximation 

('44'! rs ~ _L ~ D~ 1 -dz{S{y-h)) 1 d;,{S{y + h) - d;,{S{y)) '\ 

^ 2\dy ^ dz J ^ 2\dy J ^ 2\2 Sy{y) - Sy{y - h) ^2 Sy{y + h) - Sy{y) J 

can be used to approximate the strain. The reported value of €yz at each time step is then the average of 
the strains calculated over each tuple of tracers. Since the entire upper plate moves at a single velocity at 
any given time, ddyjdz is negligible in this case, whereas in general, this term is required to compute the 
shear strain. 

4.1.2. Computation of Stress Induced by Fluid Motion. Erom Newton’s viscosity law the 
component of stress can be found as 


(4.5) 


4 . = c(x) (^ ^ 


Since the velocity field is already known from solving the Navier-Stokes equations at each time step, the 
relevant derivatives can be approximated by finite difference approximations. As with the strain calculation, 
the second term, duyjdz^ is zero since the y velocity on the entire top plate of the rheometer is zero. The 
reported value of at each time step is then found by spatially averaging over the top 2.5 ym of the 
domain. This is done instead of just averaging over the very top of the domain in case there are numerical 
boundary layers in the fluid flow field near the boundary. Boundary layers of thickness 0{^/yAt/p) are 
known to sometimes arise in projection method based fluid solvers [31101 .To mitigate this problem, we use 
boundary conditions that do not cause this issue in the constant density and viscosity case. 

4.1.3. Computation of Stress Induced by the Biofilm Configuration. In order to compute the 
force exerted by the biofilm connections on the top plate, we integrate the Eulerian force field induced by 
bacteria adhered to the top plate. To determine if a bacteria is adhered, we choose a distance, 7 = 0.4 pm 
from the top plate, and assume that each bacteria with y coordinate in the interval [yL — 7 , yh] is adhered 
to the top, and that its ^-component of velocity is fixed to be that of the upper plate. Eor these bacteria, 
any force applied on them by spring-like connections to other bacteria behaves like a force exerted by the 
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biofilm on the upper plate instead of on the bulk fluid. The sum of these forces is used to compute the stress 
induced by the spring-like connections by means of Cauchy’s traction law, 


(4.6) 


CT n = 


p 6 


The outward unit normal, n, is (0, 1, 0) in this case since the top plate is parallel to the xz plane. The force, 
is found by integrating the Eulerian force density field that would be generated by the biofilm nodes 
adhered to the top plate. Additionally, since we are interested in the applied shear stress, this can be 

pib 

found as Note that 7 was chosen arbitrarily, however we observed that with 7 = 0.7 (im the results 

were not significantly different. 


4.2. Shear Moduli, G'and G". When a nearly isotropic material is subjected to an oscillatory dis¬ 
placement field with frequency z^, we may write the strain as e(t) = where i is the imaginary unit 

and eo is the strain amplitude. For cases where the strain is primarily only shear strain equation (4.1) gives, 
cr{u) ^ Gl{u)e{u) where Gl{u) is related to the Fourier transform in time of Gi{t). In general G^(z/) is 
a complex valued function. Breaking the complex shear modulus into its real and imaginary components, 
GJ(z^) = G'(z^) + iG"{v)] and given a strain amplitude eo(z^) and stress amplitude cro(z^) (in Pascals), 


(4.7) 


/o// \ ^o(^) c/ \ 

G (z/) = — ——co^oyy), 
eo{iy) 


G"{u) 


crojiy) 


sm6{u). 


Here, 6{iy) is known as the loss angle, measured in radians at frequency u. In the literature, G'(z^) and G" {v) 
are often referred to as the storage and loss moduli. They correspond to the elastic and viscous components 
of a viscoelastic stress strain relationship. 

Taking the domain to be a rectangular solid, oriented as shown in Figure we assume that all fields 
are periodic in the x and ^ directions. We use the following boundary conditions: 


(4.8) 


dP 

dy 


= 0 

y=^,yL 


u(x, 0, z,t) = 0 


u{x,yL,z,t) = {0,0,Ub{t )). 


Along the top boundary, we set the z velocity to be 


(4.9) 


Ub{t) = eo 


(e^^^ - I)((e^^^ - I)cosz/t + 8 e 2 ^^ sinz/t) 
(I -h 


This particular function is chosen since it is continuous, at t = 0, = 0, and because it converges to within 

0.001 of eo cos ut within half an oscillation, reducing the amount of time needed to run simulations. To 
initialize the bacteria positions, we take a 9 /im x 27 /im x 9 /im subset of a 30 /im x 30 /im x 10 /im bacteria 
position data field obtained experimentally. This data is also used in the initialization of the viscosity and 
density fields present in the biofilm. We believe that setting the internal forces to zero at the start is 
reasonable since experimental results from SAR under both compression and tension yielded similar results. 
In Figure [rTj the deformation induced by an oscillatory shearing motion is depicted. 

In order to tune our model to the experimental data, we adjusted the spring constant, kij used in Hooke’s 
Law and the distance by which we allow any two biofilm nodes to be connected by at time t = 0. For a 
spring connecting two points in space, Hooke’s Law can be written as 


(4.10) 


Fy — kij Ajj(X, Xq) (Xj — Xj), 


with 


(4.11) 


Ai,(X, Xo) 


X,(t)-X,(t)||-||X,(0)-X,-(0)|| 

||X,(t)-X,-(t)|| 


Following Hammond et al. we choose each kij to be a force constant, Fj^ax divided by the initial 

separation of bacteria i and j. Since the immersed boundary method requires a Lagrangian force density, 
we then divide Fij by the Lagrangian volume element, d^. Additionally, we note that in Dan Vo et al. [9] 
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Biofilm Connectivity at t = 0 


Biofilm Connectivity at t = 0.0250 




Figure 4.1: Starting at t = 0 on the left, the images show how the biofilm is moved as the top plate 
oscillates. Dots are bacteria locations and lines indicate viscoelastic connections. In the simulations the 
domain is periodic in the x and z directions. The periodicity is not shown here since it makes it more 
difficult to visualize the effect of deformation on the biofilm. 


and Peskin [24], an identical constitutive relation is derived from the starting point of energy functionals in 
which the force density is found by taking a Frechet derivative of an energy functional. 

In Figure 4^ we depict the frequency dependence of G' and G" . From these results, it is clear that 
our model fits experimental data on G' quite well. For G" the fit is not as strong, although we still do see 
that many of the results from simulation are within the range of experimental error. We observe that in fact 
the slope of G" is steeper than the experimental measurements. Although at this time, the cause of this 
difference is unknown, it is possible that extensions such as those discussed in Section may correct this. 


4.3. Creep Compliance Measurements J(t). Creep compliance is a measure of how a material 
deforms over time in response to an applied stress. Experimentally, the compliance is measured by using 
the rheometer to apply a step change to the shear stress on the upper plate and observing the resultant 
shear strain. A step change in stress can be written as, d(t) = where do is the magnitude of the 
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Figure 4.2: A comparison between experimentally measured results for G' and G" is shown in comparison 
the simulation results. For these results, the force constant was Fmax = 1-3223 • 10“^ , the connection 
distance between bacteria was 1.62/im, no damping was used, = 250/ro, and = 1.12po- The dashed 
lines indicate the experimental error range. 
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step change and, H{t) is the Heaviside step function. For a linear isotropic material with o-yz{t) = cr(t), the 
integral in Equation ( |4.2[ ) simplifies to 

(4.12) eyz{t) = croJi(t). 


From a physical standpoint, most of the bacteria and the bulk of the fluid are only effected by the step 
change in stress after the stress propogates vertically through the biofilm. However, in the portion of the 
biofilm adjacent to the upper plate the effect of a change in stress is instantaneous. Thus, we can write a 
force balance between the forces in the biofilm, the acceleration of the top plate, and the applied force on 
the top plate. This leads to an impulse boundary condition which specifies the velocity at the top plate The 
boundary condition can be written as 


(4.13) 



y=H 


(pV) '^{ao - - cr^)A 


where and are the stress exerted by the fluid and the springs in the biofilm at the top plate, A is the 
area of the upper plate, and {pV) is the mass associated with the top plate of the rheometer. We assume 
that this mass is equivalent to the mass of the top 2 A pm of the biofilm where bacteria are adhered to the 
top plate. 

Numerically, this boundary condition can be written as 


(4.14) \y=H = \v=H + ( -— ) (cro - ct'’ - cr^). 

ypOUoL pijk J 

In ( |4.14[ ), po and po are the density and viscosity of water, and L is the characteristic length (in this case 
10 pm). In numerical experiments, rather than immediately impose a step in stress at time 0, we add a 

mollifier, on the applied stress, ao to mitigate any possible numerical instabilities associated 

with a discontinuous boundary condition. In this case, a = 200 is chosen to be large so that the applied stress 
approaches its equilibrium value within 0.1 seconds. This is reasonable because very short time compliance 
behavior is not generally experimentally measurable, and also a step in the stress may not actually occur 
instantaneously from the perspective of a very short time scale. Results from simulations using this boundary 
condition andtwo values of ao are shown in Figure |T3l 
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Compliance vs. Time 



Figure 4.3; The time dependence of J{t, ao) is shown above for ao = 0.1 Pa and cto = 1 Pa. It can be seen 
that the compliance levels out at a similar value as the experimental result, but levels out faster than in the 
experiments. For the compliance simulations, we used Fmax = 2.9091 • 10“^ instead of F^ax = 1.3223 • 10“^ 
due to some numerical stability issues. 


Presently, we do not propose mechanisms by which the spring-like connections may break and reconnect, 
although it is possible to incorporate such a model into our current framework as discussed in [15]. Thus, 
the long term behavior of Ji(t, ao) which likely depends on the gradual redistribution of the connectivity 
of the biofilm is not expected to be captured. Therefore, we do not simulate beyond 0.1 seconds. In the 
experimental results, there are some damped oscillations present in the creep compliance experiments. In 
Pavlovsky et al. (2013), these oscillations are found to be related to intertial effects from the rheometer itself 
and thus are not expected or observed in our simulations. 

4.4. Similarity of Material Properties Between Different Bacteria Position Data Sets. In 

Dzul et al. m, the spatial statistics of bacteria in a biofilm are studied. We show here that from data 
sets that have similar spatial distributions of nearest neighbor connections between bacteria, similar bulk 
property measurements are obtained. In our simulations, we take blocks of experimental data that are 18 /im 
wide, 9 yum long, and 27 /im high and compute the dynamic moduli of each block at a fixed frequency. . The 
graphs in Figure |T^ show the stress and strain of four different biofilms over one period of oscillation. 

From figure |4.4| and Table |4.1[ we see that in three of the four biofilm data sets, similar results are 
obtained.. We also note that the mean standard error in the experimental results for this particular test was 
3.9791 for G' and 0.8836 for G". We also suspect that if larger data sets were used, even better agreement 
would be seen in the computed values of G' and G". 

The importance of this section is in verifying that the properties we are validating can be considered 
as bulk properties. Since three out of the four data sets provided results within the experimental error 
deviation, we believe this to be strong evidence the properties we measure are bulk properties. 

4.5. In-Stream Tumbling of Biofilm Fragment. Bacterial structures exhibit a diverse range of 
interaction with fluid flow. One such interaction is the tumbling motion of aggregates in shear flow. To 
simulate this effect, we conduct simulations in which there are no bacteria anchored along the plates of the 
domain. Instead, an aggregate of bacteria is located near the middle of the computational domain and the 
upper and lower plates move in opposite directions as 



(4.15) 


Uz{x, Vl, z) = 10 


1 + e-* 
15 


= -Uz{x, 0, z). 














Shear Stress vs. Time Strain vs. Time 



Figure 4.4: These graphs show the stress vs. time and strain vs. time for 4 different biofilm samples. The 
samples are all IS jam x 27 jim x 9 jam size and contain approximately 2000 bacteria positions. 


Simulation 

I 

2 

3 

G'{p = 49.91) 

13.06 

9.18 

10.03 

G"{p = 49.91) 

-5.16 

-4.44 

-3.92 

= 49.91) 

0.376 

0.451 

0.373 


Table 4.1: Results for G"(z^) and 6{u) are shown at two frequencies for 3 different biofilm coordinate 

data sets with u = 49.91 radjs. These results show that the physical properties measured here do not depend 
solely on the exact microstructure of the biofilm, but on sometype of more large scale organization of the 
bacteria positions in space. 


The boundary velocities are scaled by 10“^ so as to limit the shear stress so that the aggregate is not simply 
torn apart. To ensure that the biofilm is not attached to the plates and is sufficiently far from the plate to 
induce a rotating, or tumbling motion, these simulations only include bacteria that are greater than 8.8 /rm 
from either plate at the start. With the physical parameters we use, the bacteria aggregation rotates and is 
deformed by the fluid shear forces exerted by the fluid [6]. 

In Blaser et al. |3], analytical results on the frequency at which a solid ellipsoid will rotate in shear flow 
are provided. For an ellipsoid with axis aligned with the direction of fluid motion, the rotational frequency 
is found as 


y ^ 27r(af + a|) 
aia2r 

where r is the shear rate and aiand a 2 are the principle axes of the ellipse undergoing rotation. In our sim¬ 
ulation (see Figure [T^ , we show that a bacterial aggregate approximated as a hydrodynamically equivalent 
ellipse will rotate at a frequency similar to the theoretically expected result. The rotational frequency of 
the aggregate is found by computing the average frequency of rotation of bacteria in the yz plane about 
the center of mass of the aggregate. In the simulations, we observed a frequency of approximately 0.76 
seconds for an aggregate approximated by an ellipse with major axis ai = 9.167 ym and first semimajor 
axis a 2 = 0.136, and with shear rate of 20 containing 54 bacteria. The theoretical result in this case is 
0.641 seconds. These results are shown in Table [4^ Although in our case these results are not as concrete 
of a metric as the dynamic modulus and compliance results, rotational frequency has been used for model 
validation in fields such as red blood cell modeling (see El)- 

5. Discussion and Future Directions. Although we see that the hrIBM model provides a versatile 
means of simulating biofilms and can accurately capture some of the experimentally observed behavior of 
biofilms, there is still room to extend the model to allow for more general modeling of biofilms. An important 
area of research in biofilm studies is developing an understanding of the fracture mechanics and dynamics of 
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Positions of Bacteria at t =0 


Positions of Bacteria at t =0.398 


Positions of Bacteria at t =0.798 





Figure 4.5: Biofilm aggregate suspended in shear flow rotate over time. Snap shots shown at 0, 0.4 and 0.8 
seconds into the simulation. The flattening of the ellipse in response to the shear flow can be distinguished 
between the first and third figure. Several bacteria are marked red to help show the rotation of the aggregate. 
The blue lines indicates the trajectories of the marked cells relative to the center of mass of the aggregate. 
Distances are in 10s of micrometers. 


Major axis, ai 

First minor axis, a 2 

Theoretical Period 

Observed Period 

Relative Error 

1.67fim 

I.36/im 

0.645 

0.765 

+18.8% 

A.dAidm 

2.74/im 

0.595 

0.67s 

+13.5% 


Table 4.2: Comparison of theoretical and simulated rotational results 


biofilms. One way to model fracture dynamics in our biofilm model would be the inclusion of a stochastic 
model governing the connectivity of the viscoelastic links in the biofilm. Thus, we could define a probability 
based on the stress and strain in each link and the proximity of each pair of connected bacteria to allow 
for reconfiguration of the connectivity of the biofilm over time. It is also possible to model the viscoelastic 
properties of the biofilm by adjusting the constitutive model that is used to provide the Lagrangian force 
based on the nodal configuration of the bacteria. Modeling the changing connectivity of a biofilm was 
explored in [2j @1 126] . 

Another area that could be explored is the shape of the discrete S function used to approximate each 
bacteria and its associated viscosity halo. It is possible that adjusting this function may allow for more 
accurate modeling of the mass displacement induced by the bacteria bodies in the bulk fluid. Adjustments 
to the S function may also allow for the inclusion of nonspherical bacteria into the model. It would be 
interesting to see if similar results are obtained for bacteria that are different shapes. 

One possible difficulty in adjusting the discrete 6 function is the preservation of mass in the model. In 
Equation (2.2), there is no density dependence as would normally be seen with the Navier-Stokes equations 
for a variable density system. For the original IBM, this is in fact exact as described in [24]. For the hrIBM, 
there is an error however, error term is expected to be small in our situation since p(x, t) only varies by 12% 
over the domain (density of bacteria is not highly variable), all simulations are at low Reynold’s numbers, 
and the high viscosity gradients which overlap where density gradients occur make it difficult for fluid to 
rapidly travel down a density gradient. Although we believe that this approximation is accurate in the 
simulations we conduct, this error term may increase in situations with higher Reynold’s numbers. This is 
an area that may be further developed in future works. 

Another area of potential improvement is the numerical methods. Currently the scheme is 0{At) and 
0{h). In future work, the Crank-Nicholson time-stepping scheme may be used since, at least in the constant 
viscosity and density case, it could lead to O(At^) convergence as shown in Brown et al. [5|. Another 
approach is to use a predictor-corrector type of method such as those described in m and [24|. However, 


17 















even without heterogeneous material properties, obtaining O(At^) convergence in the overall IBM is more 
complicated and also depends on properties of the discrete Dirac S function. A detailed discussion can 
be found in Liu and Mori [niiiHi. The development of efficient O(At^) IBM schemes with heterogeneous 
material properties is still an area of active research. 


6. Conclusions. Based on the experimental results shown above and in m, we show that our biofilm 
model, Equations (2.1)-(2.8), can be used to determine bulk material properties of bacterial biofilms. In 
particular, we show that the model yields close agreement with experimental results from [22] in which the 
bacterium S. epidermidis was grown in a bioreactor and characterized using a parallel plate rheometer. We 
also show that suspended aggregates of bacteria in shear flow rotate with a similar period as a hydrodynam- 
ically equivalent ellipse. Another development is the uniformity of bulk material properties over different 
experimental data sets that possess similar spatial statistics. An important step in obtaining these results 
was the computation of bulk material properties from simulations. To our knowledge, our model is the 
first that can compute bulk material properties of biofilms based on direct simulation of both microscale 
connectivity of the biofilm and the heterogeneous rheology of the ECM. 

We also acknowledge a number of new research directions and extensions that can be done to improve 
results and also to allow for more flexible modeling of biofilms in different scenarios than what we have 
considered here. 
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